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Abstract 



The interaction potential energy surface (PES) of He-H2 is of great im- 
portance for quantum chemistry, as the simplest test case for interactions 
between a molecule and a closed-shell atom. It is also required for a detailed 
understanding of certain astrophysical processes, namely collisional excitation 
and dissociation of H2 in molecular clouds, at densities too low to be acces- 
sible experimentally. A new set of 23 703 ab initio energies was computed, 
for He-H2 geometries where the interaction energy was expected to be non- 
negligible. These have an estimated rms "random" error of ~ 0.2 millihartree 
and a systematic error of ~ 0.6 millihartree (0.4 kcal/mol). A new analytic 
He-H2 PES, with 112 parameters, was fitted to 20 203 of these new ab initio 
energies (and to an additional 4862 points generated at large separations). 
This yielded an improvement by better than an order of magnitude in the fit 
to the interaction region, relative to the best previous surfaces (which were ac- 
curate only for near-equilibrium H2 molecule sizes) . This new PES has an rms 
error of 0.95 millihartree (0.60 kcal/mole) relative to the the 14 585 ab initio 
energies that lie below twice the H2 dissociation energy, and 2.97 millihartree 
(1.87 kcal/mole) relative to the full set of 20 203 ab initio energies (the fitting 
procedure used a reduced weight for high energies, yielding a weighted rms 
error of 1.42 millihartree, i.e., 0.89 kcal/mole). These rms errors are compa- 
rable to the estimated error in the ab initio energies themselves; the conical 
intersection between the ground state and the first excited state is the largest 
source of error in the PES. 
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I. INTRODUCTION 



The interaction potential energy surface (PES) of He-H 2 (i.e., two hydrogen atoms plus 
one helium atom) is of great importance for quantum chemistry. He-H 2 is the simplest 
test case for theories of the interaction of a molecule with a closed-shell atom. It is highly 
desirable to have a PES for He-H2 that attains or at least approaches the "chemical accuracy" 
(of order one millihartree, i.e., better than one kcal/mol) required for reaction dynamics — cf. 
available ab initio potential energy surfaces for H 3 (molecule plus open-shell atom) 1-5 or H 4 
(two molecules) 6-9 . Previous analytic He-H 2 surfaces 10-24 were based on only a few dozen 
ab initio energies, which provided highly inadequate coverage of all but a very restricted 
region in the conformation space of He-H 2 . This paper reports a greatly-expanded set of 
ab initio energies, and a greatly improved He-H 2 PES fitted to these energies. 

The He-H 2 PES is also of particular astrophysical interest for studying He + H 2 interac- 
tions in physical conditions not accessible to experiment, namely the low densities charac- 
teristic of giant molecular clouds in the interstellar medium, where star formation occurs. 
Helium is the second-most-common collision partner for H 2 in these molecular clouds, which 
contain roughly one helium atom for every ten hydrogen atoms. Heating of these clouds 
by strong shock waves causes rotational and vibrational excitation of the H 2 molecules, and 
can lead to collision-induced dissociation of H 2 into free H atoms. The collision rates in 
molecular clouds can be so low that the (observed) forbidden (quadrupole) infrared emis- 
sion of excited H 2 molecules can induce highly non-thermal distributions over the internal 
states of H 2 (see, e.g., Refs. 25, 26, and 27). Because the mean free paths of molecules are 
thousands of kilometers at such low densities, these processes will remain inaccessible to 
laboratory experiment. Computer simulation is therefore a sine qua non in the study of the 
physics and chemistry of star forming regions. 

A. Notation 

Atomic units are used in this paper unless otherwise specified, i.e., distances are in 
bohrs (a ) and energies are in hartrees (Eh), millihartrees (mEh), or microhartrees (fiEh). 
Recall that 1 a = 0.529177 A, while 1 mE h = 0.0272114 eV = 0.62751 kcal/mol. 

Unless otherwise stated, energy values are measured relative to the energy of two isolated 
hydrogen atoms plus an isolated helium atom. Thus an energy E = 0.0 Eh corresponds to 
the H 2 dissociation energy, while an isolated ground-state H 2 molecule (plus an isolated 
helium atom) lies at E = —164.6 mEh — —103.3 kcal/mol. The minimum of the H 2 
potential energy curve of Schwenke 28 is E = —174.496 mE h = —109.50 kcal/mol; the lowest 
possible energy on the He-H 2 PES lies about 0.05 mE h (0.03 kcal/mol) below this, due to 
the He + H 2 van der Waals well (discussed in § II A 4). 

Figure 1 illustrates the notation used in this paper for the interatomic distances and 
angles. The (r, R, 7) coordinate system is the natural one to use when considering low- 
energy interactions of a He atom with a near-equilibrium H 2 molecule, e.g., the van der 
Waals well. It is less suitable for higher-energy interactions, where the He atom may be 
relatively close to one of the H atoms; He-H 2 surfaces that attempt to fit this interaction 
region typically use the (r, Ra, Rb) coordinate system. 



2 



B. Previous surfaces 

Many of the early He-H 2 surfaces were intended only to describe the He + H 2 van der 
Waals potential well and the outer (low-energy) part of the He-H 2 repulsive wa ll 10_16 ' 18 ' 19 ' 23 ) 
and were valid only for equilibrium H 2 molecules. However, several analytic He-H 2 surfaces 
were intended to be valid in the higher-energy interaction region. 

In 1974, Wilson, Kapral, & Burns 17 presented a more general He-H 2 PES, albeit with 
an extremely simple functional form (comprising the three pair-wise potentials plus an in- 
teraction term with only four free parameters). Shortly thereafter, Dove & Raynor 20 used a 
modified form of this surface (with somewhat improved pair-wise potentials) to investigate 
the collisional dissociation of H 2 by He. Both of these surfaces 17 ' 20 contain an unphysical 
"hole" at relatively short distances where the PES dropped abruptly to large negative val- 
ues. For these surfaces, the lowest point on the rim of this hole lies at 53 mEh (33 kcal/mol) 
and 87 mE h (55 kcal/mol), respectively (where zero energy corresponds to separated atoms: 
see § I A above); i.e., unphysical effects due to this hole could be encountered at less than 
twice the dissociation energy relative to equilibrium H 2 + He. In addition, neither of these 
surfaces 17 ' 20 included any He + H 2 van der Waals well. 

The 1982 surface of Russek & Garcia 21 was intended to fit only the repulsive wall; their 
ab initio energies had been calculated for relatively small He — H 2 separations for near- 
equilibrium H 2 molecules (of size r = 1.2, 1.4, and 1.6 a ). 

The 1985 rigid-rotor He-H 2 PES of Schaefer & Kohler 22 was fitted to the relatively 
accurate ab initio energies of Meyer, Hariharan, & Kutzelnigg 29 ; these had been calculated 
for relatively large He — H 2 separations for near-equilibrium H 2 molecules (r = 1.28, 1.449, 
and 1.618 a ), supplemented by less-accurate calculations over a wider range (0.9 a < r < 
2.0 a ). The PES fitted the van der Waals well and the He + H 2 repulsive wall, but was 
defined only for a limited range of H 2 molecule sizes (differing by no more than 40% from 
equilibrium H 2 ), and thus could not be used for cases where the H 2 molecule might be 
dissociated or highly excited. 

In 1994, Muchnick & Russek 24 presented a general He-H 2 PES, with 19 fitted param- 
eters. This surface was fitted to a combination of the Meyer, Hariharan, & Kutzelnigg 29 
ab initio energies (the less-comprehensive ab initio energies of Senff & Burton 30 were used 
as comparison values but not actually fitted), and the somewhat less accurate ab initio en- 
ergies of Russek & Garcia 21 . This Muchnick & Russek He-H 2 PES was a great improvement 
over previous surfaces: it was designed to be accurate in the van der Waals well and in the 
repulsive wall, and to behave reasonably in regions not constrained by ab initio data. 

Recently, Tao 31 presented improved ab initio He-H 2 energies in the van der Waals well 
and repulsive wall, for near-equilibrium H 2 molecules (r = 1.28, 1.449, and 1.618 a ). These 
ab initio energies should be significantly more accurate than the earlier energies of Meyer, 
Hariharan, & Kutzelnigg 29 or Senff & Burton 30 (and in fact show a van der Waals well that 
is several percent deeper than indicated by the earlier energies) . However, Tao 31 did not fit 
an analytic surface to the ab initio energies. 
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C. Goals of the present work 



An analytic general He-H 2 PES was desired that would not only accurately represent 
the van der Waals well, but also fit the interaction region at the "chemical accuracy" level 
required for reaction dynamics (i.e., of order one millihartree) . In order to do this, 20 203 
new ab initio energies were computed in the interaction region. Both the ground state 
energy and the first few excited state energies were computed, and the conical intersection 
of the ground state with the first excited state was mapped out approximately. A new 
analytic He-H 2 PES was fitted to these 20 203 ground state ab initio energies (and to an 
additional 4862 points generated to constrain the fit at large separations), yielding an order- 
of-magnitude improvement over previous He-H 2 surfaces. These ab initio energies and a 
Fortran program for this analytic PES (including analytic first derivatives) are available 
from EPAPS 32 or from the authors 33 . 

II. METHODS 
A. The grid of conformations to be fitted 

1. The main grid 

The main set of He-H 2 ab initio energies was computed for a grid of 16 703 conformations, 
defined in terms of the H 2 molecule size r, the distance Ra of the He atom from one of the 
H atoms, the angle 9 between these two, and in some cases the distance z% of the He atom 
from the midplane of the H 2 (see Fig. 1). The set of Ra grid values comprised {0.6, 0.7, 
0.85, 1.0, 1.1, 1.2, 1.3, 1.4, 1.525, 1.65, 1.775, 1.9, 2.05, 2.2, 2.35, 2.55, 2.75, 2.95, 3.2, 3.45, 
3.7, 4.0, 4.3, 4.6, 4.95, 5.3, and 5.7 a }. The H 2 -size grid values r were taken from the same 
set as Ra, supplemented by a set of extra values at short distances {0.5, 0.55, 0.65, 0.775, 
and 0.925 a }, and a set at long distances {6.1, 6.55, 7.0, 7.45, 7.9, 8.35, 8.8, 9.25, 9.7, 10.15, 
and 10.6 a } — Figure 2 shows a typical grid example (a slice of the grid, at r = 2.35 a ). 
For large H — H separations (namely, r > 6.1 a ), only conformations having Rb < 5.9 a 
were computed (the He atom lying more or less between the two H atoms) — in other words, 
the shortest two of the three interatomic distances were never allowed to exceed 5.9 a . An 
example of this is shown in Figure 3 (a grid-slice at r = 9.25 a ). 

The angle 9 was in general taken at intervals A9 = 7.5° (in 0° < 9 < 180°); for cases 
with Ra = 5.7 a or with r = 0.5, 0.55, or 0.65 a , intervals A9 = 15° were used instead, 
and cases with 9 > 90° were treated specially. Where the ray 9 = 120° reached the y- 
axis of Figure 1 (i.e., the He atom equidistant from the two H atoms), it was "bent" to 
follow the y-axis (since the grid needed to cover only the first quadrant of the y-z plane). 
Similarly, the rays 9 = 102.5°, 105°, and 97.5° were "bent" to follow the lines Z3 = r/8, r/4, 
and 3r/8, respectively; but in these cases, the rays were terminated when they reached the 
point where the usual angular separation RaA9 reached twice the separation between the 
"bent" rays (so as to avoid excessively closely spaced angles at large Ra)- For 9 > 120° (i.e., 
max{i? J 4, Rb} < r, with the He atom lying between the H atoms), a tighter ^-interval was 
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used if necessary to yield RaA9 ^ ARa; also, extra points were interspersed on the y-axis 
if necessary, so that Ay 3 ~ ARa on the y-axis. An example of this is shown in Figure 2a. 

2. Set of "random" conformations 

Another set of He-H2 ab initio energies was computed for 3500 random conformations. 
The distances were chosen in the ranges 0.6 a < r < 9 a and 0.9 a < Ra < 5.7 a using 
non-uniform probability distributions that had broad peaks near 1.5 a . The angle 9 was 
chosen using a uniform distribution, subject to the constraints that the the He atom must lie 
in the first quadrant of the y-z plane and that no atom could lie further than 5.7 a from its 
nearest neighbor. Also, "near-duplicate" conformations were not permitted (i.e., where all 
cartesian coordinates agreed within 0.01 a with the coordinates of another conformation). 
Figure 2a shows those random conformations that lie near the grid-slice r = 2.35 a . 

3. Unfitted "random" conformations 

The same probability distributions were used to obtain a second set of 3500 random 
He-H 2 conformations, for testing purposes. The resulting ab initio energies were used only 
to check the rms error at "interstitial" conformations; they were not used in the surface 
fitting process. 

4- Accurate van der Waals He+ il 2 energies 

A set of 2145 accurate energies was generated for conformations in the van der Waals well 
and outer repulsive wall, for H 2 molecules not too far from equilibrium. These energies were 
generated from the accurate ab initio energies and fitted surfaces of other authors 22 ' 24 ' 31 , as 
described below. These accurate van der Waals energies were given high weight in the fit, 
with lower weight for smaller R or for r further from equilibrium (as discussed in § II D 1). 

Unlike the main grid, these conformations were defined in terms of the distance and 
angle of the He atom relative to the H 2 molecule, i.e., in terms of {r, R, 7} (rather than 
{r,R A ,9}). The set of r values comprised {0.6, 0.7, 0.8, 0.9, 1.1, 1.28, 1.449, 1.618, 1.8, 2.0, 

2.4, 3.0, and 4.0 a }; the set of R values comprised {4.0, 4.5, 5.0, 5.25, 5.5, 5.75, 6.0, 6.25, 

6.5, 6.75, 7.0, 7.25, 7.5, 7.75, 8.0, 8.5, 9.0, 10.0, 11.0, 12.0, 15.0, and 20.0 a }. The angle 7 
was generally taken at 15° intervals (in 0° < 7 < 90°); for R = 20.0 a , and for roughly 
half of the R values at non-equilibrium r, intervals of 30° were used instead. Conformations 
generated at relatively small R were discarded if they had an interaction energy of more 
than a few mE h between the He atom and the H 2 molecule (i.e., only the van der Waals 
well and the very outermost part of the "repulsive wall" were included). An example of this 
grid is shown by the squares in Figure 2b (a grid-slice at r = 2.4 a ). Note that there is an 
overlap with the main ab initio grid. In the outer part of this overlap region, the interaction 
energy of He with the H 2 molecule is comparable to the uncertainty in the main-grid ab initio 
energies, so the fit depends largely on the generated van der Waals H 2 + He energies. In the 
inner part of the overlap region, the interaction energies are larger and the generated van 
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der Waals H 2 + He energies are given fairly low weight (see § II D 1), so both sets of energies 
contribute to the fit. 

Tao 31 used the complete fourth-order M0ller-Plesset approximation (MP4) 34 ' 35 with a 
large basis set to compute 69 accurate ab initio He-H 2 energies in the van der Waals well 
and repulsive wall (15 R values in 2 a < R < 15 a ), for three orientations (7 = 0°, 
45°, and 90°); separations R over the full range were reported only for the equilibrium H2 
molecule (r = 1.449 a ), but a few (R = 3.0, 5.0, 6.5, and 8.0 a ) were also reported at 
r = 1.28 and 1.618 a . (Note that Tao reported a maximum van der Waals well depth of 
47 \iE h for He plus equilibrium H 2 , at r = 1.449 a , R = 6.5 a , 7 = 0°.) This subset of 
ab initio energies was referred to as the "Tao vdW" points. 

For r = 1.28, 1.449, and 1.618 a , a relatively simple formula was fitted (for 7 = 0°, 
45°, and 90°) to the difference between the Tao 31 ab initio energies and the analytic surface 
of Schaefer & Kohler 22 , and used to interpolate corrected energies at R and 7 values not 
computed by Tao. These generated energies were referred to as "gen- Tao vdW" points. 

For r < 1.1 a and r > 1.8 a , energies were obtained from a version of the Muchnick 
& Russek 24 surface where the van der Waals term in their potential had been modified 
slightly to better fit the newer, more accurate Tao 31 ab initio energies (this surface is the 
one described in § II D 2). These generated energies were referred to as "gen-MR vdW" 
points. 

For r = 0.9, 1.1, 1.8, and 2.0 a , energies from the Schaefer & Kohler 22 surface were also 
used (since it might give energies of comparable accuracy in this region). These generated 
energies were referred to as "gen-SK vdW" points. 

5. Accurate "H-He" energies 

A set of 1887 energies was generated for conformations consisting of a relatively close- 
together H — He pair (of size Ra) with a very-distant H atom (a distance R' from the nearest 
point on the line segment Ra)- Since the contribution of the isolated H atom was completely 
negligible, this subset was referred to as "H-He" points. 

For Ra <4a , the H-He energies were generated from 31 ab initio He + H energies that 
we computed at R A = {0.5, 0.55, 0.6, 0.65, 0.7, 0.75, 0.8, 0.9, 1.0, 1.1, 1.2, 1.3, 1.35, 1.4, 
1.45, 1.5, 1.6, 1.7, 1.8, 1.9, 2.0, 2.2, 2.4, 2.6, 2.8, 3.0, 3.2, 3.4, 3.6, 3.8, and 4.0 a } (note 
that we had also computed 9 other ab initio energies at larger Ra values, which were not 
used). For Ra > 4 a , the H-He energies were generated from the H — He pair-wise potential 
of Dove & Raynor 20 (which used slightly-modified functions from Gengenbach, Hahn, & 
Toennies 36 ), with values chosen at R A = {4.25, 4.5, 4.75, 5.0, 5.25, 5.5, 5.75, 6.0, 6.25, 6.5, 
6.75, 7.0, 7.25, 7.5, 7.75, 8.0, 8.5, 9.0, 10.0, 11.0, 12.0, 15.0, and 20.0 a }. Distances to the 
isolated H atom comprised R' = {7.1, 7.7, 8.3, 9.3, 10.3, 11.4, 13.1, 14.8, 16.5, 19.0, 22.0, 
and 25.0 a }. Conformations were discarded whenever pair-wise H — He or H — H energies 
from the isolated H atom would be non-negligible compared to the energy at separation Ra 
of the relatively close-together H — He pair. (Although they were negligibly small for all 
retained conformations, pair-wise interaction energies of the isolated H atom with the other 
two atoms were nonetheless added to the main H-He energy, for consistency with the "H-He 
plus H" energies of § II A 6 below). For each pair {R A ,R'}, eight angular grid-points were 
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chosen, four with r = R! at ip = 90°, 120°, 150°, and 180°, and four with R B = R! at 
* = 90°, 120°, 150°, and 180° (see Fig. 1) — a "race-track" shaped grid. For R' > 2 R A , 
the two points at ip = 90° and \1/ = 90° were replaced by a single point at a distance R' 
from the midpoint of Ra- Higher weight was given to conformations with large Ra values 
(as discussed in § II D 1). 

6. Approximate "H-He + H" energies 

A set of 830 approximate energies was generated to constrain the fit in regions not covered 
by the above grids, where the two H atoms were fairly far apart (5 a < r < 19 a ), and at 
least one of the H atoms was well separated from the He atom so that the interaction energy 
of this isolated H atom would be small compared to the uncertainty in our ab initio energies 
(but might not be small compared to the van der Waals well depth). Since the contribution 
of the isolated H atom was small but not completely negligible, this subset was referred to 
as "H-He + H" points. 

The grid was similar to the main grid, but with r = {5.0, 6.55, 7.2, 7.9, 8.5, 9.25, 9.9, 
10.6, 11.4, 14.8, and 19.0 a } and R A = {0.8, 1.1, 1.4, 1.7, 2.1, 2.5, 2.9, 3.3, 3.7, 4.2, 4.7, 5.2, 
5.7, 6.2, 6.7, 7.3, 8.05, 9.0, 10.0, 12.0, 15.0, and 20.0 a }; the minimum angular spacing was 
15° in 9. Many parts of this grid were in fact more coarsely spaced than this in Ra and/or 8, 
and for r = 7.2, 8.5, 9.9, and 11.4 almost all gridpoints with 9 < 105° were omitted. In 
addition, conformations relatively close to areas covered by the other grids were omitted; 
this restriction eliminated cases where more than one of the three pair-wise potentials would 
exceed ~ 0.1 mE^. Approximate energies were thus computed by simply summing the 
three pair-wise potentials. Cases where this computed energy would be expected to be quite 
accurate (i.e., large Ra) were given increased weight (as discussed in § II Dl). 

B. Ab initio computations and analysis 

1. Computational methods and CPU-time 

The ab initio computations and analysis of errors largely follow the methods described 
in our previous papers 37 ' 38,3 ' 4,8 ' 9 . Energies were computed using the workstation version 39 of 
Buenker's MRD-CI program 40 . At the position of each of the two hydrogen atoms, we used 
the (9s3p2d)/[As3p2d\ Gaussian basis set from Boothroyd et al. 37 ' 38 (with rf-functions opti- 
mized for H 2 ). At the position of the helium atom, we used a (10s4p2rf)/[5s4p2(i] Gaussian 
basis set, with s-functions taken from a (partially-decontracted) set from Huzinga 41 , and 
the p- and ci-functions chosen based on hints from several papers 41-43 ; the basis set error 
for an isolated He atom is only 1.028 mE h with this basis set. These two basis sets are 
given in Table I; they comprise the largest combination of basis sets allowed by array sizes 
in the MRD-CI program 39 . For most conformations, molecular orbitals were obtained from 
closed shell SCF; if closed shell SCF iterations were slow to converge, open shell SCF or 
mixed (1 closed, 2 open) shell SCF was automatically used instead. A configuration selection 
threshold of T = 2 jiEh was used for most points, with T = 10 fiEh used for 3006 high- 
energy cases at short distances (Buenker's MRD-CI program automatically extrapolates to 
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zero threshold, using truncated-CI energies at T and at 2T plus a parameterization of the 
estimated effects of neglected configurations). Multiple CI roots were computed, to ensure 
that the ground state energy was always among the energies reported by the program. A 
fairly extensive set of reference configurations was used for the CI calculation: the minimum 
size of C 2 (the sum of reference configuration C 2 values) was 0.969, and the average size 
was 0.982. (Note that additional reference configurations were included automatically if the 
program estimated that they would increase C 2 by more than 0.001 individually and by 
more than 0.003 collectively). 

The above relatively small basis set and relatively large selection threshold were required 
by CPU-time considerations. So long as the ab initio energies are of accuracy comparable 
to the desired "chemical accuracy," the most important factor in obtaining a good fitted 
analytic surface is adequate coverage in the 3-dimensional conformation space of He-H2, 
which requires large numbers of ab initio energies. In total, over a year of CPU-time plus 
1/ O-time was required for the total of 23 406 ab initio computations (this includes 3203 points 
that were initially computed with only a single CI root, or with T = 10 /iE h , and which were 
subsequently re-computed more accurately). These computations were distributed among 
about a dozen different computers, mostly workstations of various types, but including a 
couple of higher-performance machines. For the faster computers the I/O-time tended to 
dominate (due to large temporary files), but there was one exception to this: over two thirds 
of the energies were computed on a single higher-end SGI with a fast disk, which allowed 
computation of a typical ab initio energy in about 10 minutes. The other computers typically 
required 30 to 60 minutes per point, with some older machines requiring up to several hours 
per point. 

2. Small corrections applied to the MRD-CI energies 

Extrapolation to zero threshold (performed automatically by the MRD-CI program) and 
the Davidson correction to full CI are both standard procedures. Although the MRD-CI 
program uses only the energy at thresholds T and 2T to extrapolate to zero threshold, it also 
calculates the energy at 3T and 4T; we therefore used these latter to make a very small (rms 
size 0.05 mE h ) "curvature correction" to the extrapolation (as Boothroyd et a/. 37 ' 38 ' 8 ' 9 did 
for H 4 — but unlike H 4 , this was the only modification of the extrapolation). The average 
size of the extrapolation was 1.73 mEh for the threshold T = 2 fiEh, and about 3 times as 
large for T = 10 fiEh- For the 17025 cases where energies were obtained with both T = 10 
and 2 fj,Eh, the rms difference between extrapolated energies was 0.17 mEh- This should 
be an overestimate of the extrapolation error of either, suggesting that the rms error in the 
extrapolation to zero threshold is about 0.15 mE h for T = 10 fiE h (used for 3006 high-energy 
points) and about 0.1 mEh for T = 2 fiEh (used for the other 17197 points). This error 
estimate is also supported by the fact that, for 996 points where a single-root T = 0.4 fiEh 
case was also computed, the rms difference was 0.12 mE h relative to the 5-root T = 2 fiE h 
cases. 

As in Boothroyd et a/. 37,38,8 ' 9 , we made a parameterized Davidson-type correction to 
full CI, namely, 

A full _ci = \% C C F) ^E DC = \% c c p \l-C 2 R )AE sd /C 2 R , (1) 
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where AE s d is the contribution of single and double excitations to the energy, and AEdc 
is the standard Davidson correction. For single-reference closed-shell configurations with at 
least four electrons, if the effects of quadruple excitations are assumed to be much larger 
than triple excitations, then the above formula for the Davidson correction (with an expected 
value of \dc ~ 0.5 for closed-shell He-H 2 ) can be derived theoretically 44-47 . We not only 
used multiple references, we also used open and mixed shell SCF cases; we therefore obtained 
values for the parameter Xoc by minimizing the average differences between MRD-CI energy 
values computed for the same point with different reference sets (which thus had different 
sizes for the Davidson correction). We used A^' = 0.52, 0.16, and 0.32 for closed, open, 
and mixed SCF cases, respectively; the average size of AE DC was 1.31, 1.90, and 0.60 mE h , 
respectively, yielding full-CI corrections averaging 0.68, 0.30, and 0.19 mi?/,, respectively. 
For the 12 077 cases where the program obtained improved reference sets, and thus two 
different-sized Davidson corrections were available, the rms difference between the resulting 
energies was 0.16 mE h . This may be a slight underestimate of the energy uncertainty 
resulting from the Davidson correction, since on average the initial reference set yielded a 
Davidson correction only 1.75 times as large as the improved reference set, and we have 
no information from cases where only one value of the Davidson correction was available. 
There are suggestions of a systematic uncertainty of order 0.1 in the preferred value of A^*^, 
leading to a systematic uncertainty in the energy correction of order 0.13 mE h . We estimate 
the total uncertainty in the Davidson-type correction to full-CI to be about 0.3 mE h . 

The final energy value Ea na \ was obtained from the estimated full-CI energy using a 
London-type basis correction similar to that for H4 described in Boothroyd et a/. 37 ' 38 ' 8 ' 9 . 
The H4 London-type basis correction is obtained by taking the difference between the H4 
London energy with the accurate H 2 singlet curve and the H4 London energy with the H 2 
singlet curve computed using the finite basis set. For He-H 2 , a London- type basis correction 
was obtained by taking the H 4 basis correction, letting one distance go to zero (to yield the 
He atom), and replacing the contribution from the zero distance by the isolated He atom 
basis correction of 1.028 mEh- (After doing the algebra, this turns out to mean taking the 
H 2 singlet basis correction for the H 2 molecule, adding the isolated He atom basis correction, 
and adding half of the H 2 singlet basis correction for each of the H — He distances.) This 
yields the correct basis correction in the limit of a separated He atom plus H 2 molecule. The 
average size of this correction was 2.32 mEh (or 1.29 mEh if one excludes the He atom basis 
correction). Boothroyd et a/. 8,9 concluded in the case of H4 that the systematic error in 
the basis correction might be as much as 1 mEh- For He-H 2 , the error should be somewhat 
smaller; one might estimate that the error could be comparable to the non-H 2 , non-He part 
of the basis correction, i.e., about 0.6 mE h - Comparison of a number of energies that were 
also computed using smaller basis sets suggests that this error estimate is reasonable, as 
does comparison with energies computed by Tao 31 using a larger basis set. The error in the 
basis correction will be larger for geometries where the He atom is close to one (or both) 
of the H atoms, and smaller for geometries where the He atom is relatively far from the 
H atoms. 
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3. Best MRD-CI energies and estimated errors 



Where more than one energy value had been computed for a given conformation, a 
weighted average of the "best" values was obtained for use in the fitting procedure. (In 
general, multiple-root cases were preferred to single-root cases, closed shell cases to open 
and mixed shell cases, lower thresholds to higher ones, and smaller Davidson corrections 
to larger ones.) The extrapolation to zero threshold and the Davidson- type correction to 
full CI appear to yield largely "random" errors; the systematic components of these, though 
not completely negligible, are much smaller than the systematic error in the basis correction. 
Thus the total uncertainty in the ab initio energies probably comprises a "random" error of 
about 0.2 mEh and a systematic error of order 0.6 mEh- The ab initio energies are available 
from EPAPS 32 or from the authors 33 . 



C. Functional Representation of the He-H2 Surface 

This section lays out the equations which underlie our fits to the He-H 2 ab initio data. 
Since the Muchnick & Russek 24 surface was the best of the previous surfaces, we based our 
analytic formulae on theirs, but added terms to get more flexibility. The surface is defined 
in terms of the three interatomic distances 1Z = {r, Ra, Rb}, where r is the H — H distance 
and Ra and Rb are the two He — H distances (see Fig. 1). As in Muchnick & Russek 24 , new 
(basically spheroidal) variables are defined as 

- Ra + Rb Ra — Rb /ri x 

2~ ' V = —^- (2) 

Note that R — > R and rj — > cos 7 for R ^> r, i.e., the usual Legendre function expansion 
terms Pf,(cos7) / R n in the van der Waals well can be replaced by PL(rj)/R n . However, the 
functional form of the fitted surface does not explicitly contain either R or 7. 

1. Functional form of our fitted surfaces 
The full He-H 2 potential energy surface is given by 

VkeH 2 (^) = Vk a (r)+Vi nt (^), (3) 

where 

v int (K) = Vc(n) + v H {H) + v cb (n) + v d {1i) + v D (n) . (4) 

The Vh 2 term is the isolated-H 2 -molecule potential (which is of course a function only of the 
distance r between the two H atoms); for this we used the accurate H2-molecule potential 
of Schwenke 28 . The interaction energy Vi n t between the H 2 "molecule" and the He atom 
was fitted by several terms, functions of all three interatomic distances (all terms except Vp 
being generalizations of terms used by Muchnick & Russek 24 ). The basic forms of these 
terms correspond to Coulomb, Pauli, and dispersion terms, as discussed by Muchnick & 
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Russek 24 . However, the added flexibility (provided by the summations in the equations 
below) largely blurs the distinctions between these types of terms in our fitted surfaces. 
The Vq term has the form of a Coulomb term (omitting the H — H contribution): 



V C (K) = A C {K) 



1 + 



A C {K) 



F c (Ra) + F C (R B ) 



(5) 



where 



F c (Rv) 



R v 



^ Ft 3 

.7 = 1 nv 



for v = A,B 



(6) 



and 



A c {n) 



m=l L=o i=Q 



M c 



^ a c {m,L,i) 



L=0 
L even 



R i 



(7) 



note that the lower limit m = 1 on the summation (rather than m — 0) was chosen in order 
that Ac(TZ) — > 2 for r — > as well as for r — > oo (as in the similar term of Muchnick & 
Russek 24 ). Parameters to be fitted include \c, flc, an d all of the ac(m,L,i) and bc(j)- 
Our final adopted surface had {M C) L C) Ic, Jc} — {2, 2, 1, 5}, but we tested fits with ranges 
as high as {3,4,2,5}. Note that Muchnick & Russek 24 used this V c term with only the 
parameters \c, /3c, and ac(l,0, 0). 

The Vh term looks like a Pauli-type interaction between the He atom and each of the 
H atoms: 



v H {H) = a h 



r /3 H r ^ r m P^(t7)^ a H{rn,i^,i) 

m=0 L=o i=Q 

L even 



R i 



F H (R A ) + F H (R B ) 



(8) 



where a# (0, 0, 0) = 1 is a constant (in order that the lowest-order form of Vh vanish at 
r = 0, for consistency with the similar term of Muchnick & Russek 24 ), and 



F H {Rv) 



7=1 ^ 



for v = A,B . 



(9) 



Parameters to be fitted include Ah, Xh, (3h, and all of the bc(j) and ac(m,L,i) except for 
a_ff(0, 0, 0). Our final adopted surface had {M H , L H , Ih, Jh} = {2, 4, 2, 5}, but we tested fits 
with ranges as high as {3, 6, 3, 5}. Note that Muchnick & Russek 24 used this Vh term with 
only the parameters Ah, Xh, and (3h- 

The V c b term looks like a Pauli-type interaction between the He atom and that part of 
the electronic distribution drawn in to form the H — H covalent bond: 



V cb {H) = e-Aar-A^r)* J2r m J2 Pl(v) £ a cb (m,L,t)R , 

m=0 £=0 ^=/ cfe 

L even 



(10) 



where 
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A cb (r) = e cb + ( cb e~ 5 ^ r . (11) 

Parameters to be fitted include (3 cb , e cb , ( c b, &cb, and all of the a cb (m, L, i). Our final adopted 
surface had {M cb , L cb) I cb) J cb } = {2,4,-1,2}, but we tested fits with ranges as high as 
{3,6, —2,2}. Note that Muchnick & Russek 24 used this V cb term with only the parameters 
f3 c b, £cb, (cb, S cb , a cb (0,0,0), and a cfe (0,2,0). 

The V d term is a dispersion (long-range) interaction between the He atom and the H 2 
molecule: 



V d (n) = -f Ad A d (r) 



+ 



M d 



l + fa d a d (r)P 2 (rj)+J2e £{(i ' L)r PL(v) £ a d (6,L,m) 



L=4 m=l 

L oven 



R Q + Ci 6 

C 3 [l + f ad a d (r)P 2 
R 8 + C 2 8 



+ £ ^ T lr (r))r £ fWPLlri) £ a d (n,L,m)r* 

n=g JrC + [^ d {n)\ L=0 m =min{L,l} 

L even 



, (12) 



where 



A d (r) = (4.790521 r 2 -2.494444 r + 3.008850) e ~°- 589ir (13) 

and 

a d (r) = (0.238248 r 2 - 0.126241 r) e -°- 8075r (14) 

(where r is in a ). Note that the values of e(n, L) were close to —1.0 in all fits, so that the 
term containing C 2 and C3 (retained for consistency with the form of Muchnick & Russek 24 ) 
has a significantly different behavior from the n = 8 term in the last of the summations. 
Note also that only the L = terms are allowed to have m — 0, since the angular variation 
in V d must vanish for r — > 0. For cases (such as the final fitted surface) with an upper limit 
L d < 4, the last term in the large square brackets is omitted (since L = 4 is the lowest L 
value allowed in that term). 

The 7 parameters quoted in the terms A d (r) and a d (r) above were obtained by fitting 
(with an accuracy of a few percent) to the accurate values calculated by Thakkar et al. 23 as 
a function of r for these asymptotic terms, as shown in Figure 4. The remaining parameters 
to be fitted include fA d , fad, Ci, C2, C3, and all of the C d (n), e(n, L), and a d (n, L, m). Since 
the parameters f'A d and f ad multiply a version of the asymptotic long-range form (although 
"softened" at small R by the parameter Ci), these parameters might be expected to remain 
within a few percent of unity in a fit; this did indeed prove to be the case, particularly for f^ d 
(note that a d is much smaller than A d , and thus f ad is determined with less precision by 
the fitting process). 

Our final adopted surface had {N d ,L d ,M d } = {8,2,1}, but we tested fits with ranges 
as high as {10,4,2}. Note that Muchnick & Russek 24 used this V d term with only the 
parameters C±, C 2 , and C3; they also replaced P 2 {rf) with rf, and set A d (r) = 0.81 + 1.92 r 
and a d (r) = 0.112 + 0.06 r (although they noted that these latter two functions really 
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should have used a decreasing exponential in r). Their versions of Ad(r) and a d (r) are also 
illustrated in Figure 4, although their a d {r) is not strictly comparable, since it multiplies 
rf rather than i^^)- ( in spite of the fact that their A d and increase without limit with 
increasing r, the V d term used by Muchnick & Russek 24 remains bounded, since by definition 
R > r/2, and thus their form of V d is proportional to r~ 4 for r 3> Ci.) 

The Vd term is a long-range interaction between the He atom and separated H atoms: 

V D {H) = S D (r) [F D (R A ) + F D (R B )] , (15) 

where 

MR*) = E p n a in\^n for u = A,B (16) 



n=6 
n even 



R„ n +[C D (n)f 



Mr) = 1 - ~ m ' , (17) 



and Sr>(r) is a switch function to turn Vd on with increasing r: 

£?p(r 
{6^ + [^(r)] 2 } 

with 

S D (r) = 4.790521 (r 2 + 6) e -°- 589ir (18) 

(where r is in a ). The softening parameter of 6 in the switch function Sn(r) above (which 
regulates the switchover at relatively small r) may be considered to be a surface parameter, 
but the other constants in the switch function were taken from Ad(r) in Equation (13), so 
that the Vd term would be turned on with increasing r as the Vd term was being turned 
off by A d (r). Other parameters to be fitted include all of the C D (n) and a D (n). Our final 
adopted surface had N D = 8, but we also tested fits with N D = 10. Note that Muchnick & 
Russek 24 did not include any such Vd term. 



2. Other types of terms tested 

Various index ranges were tried for the summations in the above equations. In addition, 
it was tested whether a better fit could be obtained if the terms V c , V H , and V ch were 
multiplied by an exponential cutoff in the overall size of the He-H 2 conformation, namely, 
by a factor 

e~ ppa , where p = (r 2 + R A 2 + R B 2 ) 1/2 (19) 

(with p and a being parameters that could be fitted). It turned out that such a cutoff did 
not improve the fit, and it was discarded. 

It was tested whether adding a many body expansion term V M to the He-H 2 PES would 
improve the fit: 
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Km min{fc,MM— 1— fe} min{/M,A^M— j— k} 

M#) = E E [paW+PaW] E a M {i,j,k)p r \ (20) 

fe=l j=0 i=l 



where 



p r = re 



-P M r 



and 




for v — A, B , 



(21) 



with parameters to be fitted including j3 M and all of the a,M(i,j,k). This Vm term could 
also be multiplied by an exponential cutoff as given by Equation 19. Various orders M M 
and individual index limits Km and 1m were tested, but Vm always turned out to be of very 
little help in improving the fit (compared to the other terms of § II C 1), so it was discarded. 



This section outlines in general the steps followed in developing the terms and optimizing 
the parameters in our analytical representation of the He-H 2 surface. The details of the 
equations are given in § II C 1 above. 



In this section, the weights referred to are those that were used to multiply the deviations 
between the fit and the data, before squaring these weighted deviations to calculate the rms 
of the fit. Note that, frequently, the "weight" in a fit is defined as that applied to the square 
of the (unweighted) deviation, i.e., with this latter definition the weight would be the square 
of the values reported below. 

Several different criteria were used to determine weight values. If more than one applied, 
the final combined weight used was the product of the individual weight factors from the 
following separate criteria. 

High-energy portions of the surface are less likely to be accessed in collisions of a hydrogen 
molecule with a helium atom; also, higher ab initio energies have larger uncertainties than 
lower energies. Thus all points with high energy E were given reduced weight, namely, 



i.e., a weight inversely proportional to the energy E for cases with energies more than 
about 2.2 times the H 2 dissociation energy above that of an equilibrium H 2 molecule plus a 
separated He atom. 

For non-compact conformations where the He atom was relatively far from either H atom, 
the (absolute) size of the uncertainties in our ab initio energies was expected to be smaller 
than for more compact conformations (especially for the error in the basis correction: see 
§ II B 2). Thus our ab initio energies (described in § II A 1 and § II A 2) were given an 
additional weight factor 



D. Our Approach to Fitting 



1. Weights applied to fitted energies 




, E < 0.2 E h 
, E>0.2 E h 



(22) 
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1 

(R m — 2) 
2 



Rm ^ 3 Oq 

3 a < -R m < 4 a 



where _R m = minjiiU, Rb} (23) 



(i.e., a weight increasing linearly from 1 to 2 as the shortest He — H distance R m is increased 
from 3 to 4 a ). 

The van der Waals He + H2 points described in § II A 4 are expected to have smaller (ab- 
solute) uncertainties as the distance R between the H 2 molecule and the He atom increases; 
they should also be most accurate for near-equilibrium H 2 molecule sizes (i.e., r ~ 1.4 a ). 
For these "vdW" points, two weight factors w R and w r were applied: 



w R (R) 



' 0.3 


, R < 3 a 


1 


, R = 4 a 


10 


, R — 5 a 


300 


, R — 6 a 


1200 


, R = 6.5 a 


k 3000 


, i? > 10 a 



w r (r) 



0.1 
0.3 
1 
1 

0.3 
0.1 
0.01 



r < 0.6 a 
r = 0.9 a 
r = 1.2 a 
r = 1.62 a 
r = 2.0 a 
r = 2.4 a 
r > 4.0 a 



(24) 



(where wr and w r were interpolated linearly in between the distance values specified above). 
The "Tao vdW" points (ab initio energies of Tao 31 ) were given an additional weight factor 
wt = 3, and the "gen- Tao vdW" points were given an additional weight factor w 9 t = 2. 
Note that quite a large number of fits were performed where wr was either increased or 
decreased relative to Equation (24) by about a factor of 3 for R > 5 a , but these fits 
were discarded — the weight wr of Equation (24) appeared to work the best, neither over- 
emphasizing nor under-emphasizing the van der Waals well energies relative to the ab initio 
energies. 

The "H-He" points described in § II A 5 should be as accurate as any of our other ab initio 
energies for R^ <5a (where Ra is the shorter H — He distance for these points: Ra < Rb)', 
they should be rather more accurate in the region of the H — He van der Waals well and 
further out. Thus, a weight factor WHHe(RA) = ^s^{wr(Ra), 1-0} was applied to them, 
where wr has the form given in Equation (24). 

The "H-He + H" points described in § II A 6 are expected to be somewhat less accurate; a 
weight factor WHHe+H(RA) = max{0.01 Wr(Ra), w„ c (-Ra)} was applied to them, where w nc 
has the form given in Equation (23). 



2. Modified Muchnick & Russek surface 

The first step was to get a very slightly modified version of the Muchnick & Russek 24 
surface, which we refer to as the "modMR" surface. This differed from the surface of 
Muchnick & Russek 24 in that we used the lowest-order version of our functional form 
for Vd, including use of our Ad(r) and ad{r) from Equations 13 and 14, and the fact 
that ttrf(r) multiplied P 2 (^) rather than rf . We also added the lowest order Vd term 
(i.e., N D = 6) from Equation (15). Since we used our already- fitted versions (8 param- 
eters) of A d (r), a d (r) and S D (r) from Equations (13), (14), and (17), this initial "modMR" 
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surface had 19 other basic parameters that could be fitted. There were 9 short-range 
non-linear parameters {Ac, (5c, 0, 0), Xh, (5h, Pa, £ cb, Ccb, Scb}, 3 short-range linear pa- 
rameters {A H , a c b(0, 0, 0), a c b(0, 2, 0)}, 3 long-range non-linear parameters {C\, C 2 , (7^(6)}, 
and 4 long-range linear parameters {fAd, fad, C3, d(6)}. For this first surface, we obtained 
od(6) from Gengenbach, Hahn, & Toennies 36 , and set Cd(6) = 6.7876662396 a (actually, 
2 1 / 6 x 3.2 A — a reasonable softening value to keep the Vd term from contributing signifi- 
cantly at relatively compact geometries). We set the other parameters equal to the values 
used by Muchnick & Russek 24 , and fitted only the parameters fAd and f a d, using only the 
"Tao vdW" and "gen-Tao vdW" points described in § II A 4. This procedure yielded values 
of fAd = 0.98234 and f ad = 1.261, reasonably close to unity (i.e., at large H 2 — He sepa- 
rations R, this surface still had a form very close to that calculated by Thakkar et al. 23 ). 
The surface was almost identical to that of Muchnick & Russek 24 , but was a significantly 
better fit to the more recent (and more accurate) van der Waals ab initio energies of Tao 31 
(unweighted rms error of 1 uE h for these points, rather than 3 fxE h ). This "modMR" surface 
was the one used to generate the "gen-MR vdW" points described in § II A 4. 

3. Optimization of fitted parameters 

Several hundred fits were subsequently performed to the full set of 20 203 ab initio ener- 
gies and 4862 generated points, with various ranges for the indices in the summations in the 
surface equations of § II C 1; the number of parameters actually fitted varied from a couple of 
dozen to about 150. Fitting was performed using the NAG Fortran library non-linear fitting 
routine E04FDF to minimize the weighted rms error with respect to the 20 203 ab initio 
plus 4862 generated energies, using the weights of § II D 1 (the covariance matrix was then 
produced by the NAG routine E04YCF) — a typical fit required a few hours on a higher-end 
computer. The first fit described in the previous paragraph supplied initial values for the 
basic parameters in the next fits, often with values shifted randomly; the added parameters 
in the summations were either initialized to zero or to random values. A reasonably good 
fit was frequently used as the basis for subsequent larger fits (sometimes with random shifts 
applied to the parameters before fitting). There were some exceptions to this: the 3 basic 
long-range non-linear parameters {C\, C2, C_d(6)} were not re-fitted, and other long-range 
non-linear parameters were generally not fitted, but just given a reasonable initial value. 
Also, 3 of the basic long-range linear parameters, namely, {fAd, fad, C3}, were not re-fitted 
until near the end of the fitting process. Note that in our final adopted fit the values of 
fAd ~ 0.9552 and f a( i ~ 1.149 were near unity; since these parameters multiply forms Az( r ) 
and adij-) that had been fitted with an accuracy of a few percent to the values of Thakkar 
et al. 23 (as discussed in § II CI), this final adopted fit yielded an asymptotic interaction 
between the H 2 molecule and the He atom very close to that calculated by Thakkar et al. 23 

4- Selection of a good fitted PES 

For a quick comparison of the quality of all the different fits, we considered the full 
weighted rms minimized by the fitting program with respect to all 24 804 fitted points, and 
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also the energy-weighted rms for 9 subsets of these points, where only the weight we from 
Equation (22) had been applied. These 9 subsets comprised: 

1. all 20 203 ab initio energies, 

2. the 14 585 ab initio energies with E < 0.174 E h (i.e., below about twice the H 2 
dissociation energy), 

3. the 7177 ab initio energies below the H 2 dissociation energy (E < 0), 

4. 931 ab initio energies with near-equilibrium H 2 (1.2 a < r < 1.62 a ) and a not-too- 
distant He (2 a < R < 4 a ), 

5. 468 fairly-close "H-He" points (with 2 a < R A < 4 a ), 

6. 143 van der Waals "H-He" points (with 6 a < R A < 10 a ), 

7. 231 "Tao vdw" and "gen-Tao vdw" H 2 +He van der Waals points with near-equilibrium 
H 2 (i.e., 1.28 a <r < 1.62 a , and 6 a < R < 10 a ), 

8. 434 "gen-MR vdW" and "gen-SK vdW" van der Waals points with small-H 2 (r < 1.2 a 
and 6 a < R < 10 a ), 

9. 714 "gen-MR vdW," "gen-SK vdW," and "H-He + H" van der Waals points with 
large-H 2 (r > 1.62 a and 6 a < R < 10 a () ). 

For a few cases where the above rms values suggested that the surface was among the best of 
the fits, a collection of weighted, energy-weighted, and unweighted rms values was considered 
for a more comprehensive list of about 100 subsets. 

The fitting program produced estimated uncertainties of the fitted parameters (from the 
NAG covariance routine E04YCF). When an enlarged fit did not significantly reduce the 
rms error and yielded new parameters whose values were not significant (i.e., comparable 
to or smaller than the estimated errors in these parameters), this suggested that the point 
of diminishing returns had been reached — larger fits would presumably yield spurious 
"wiggles" rather than an improved fit. This criterion was used to obtain a near-optimum set 
of index ranges for the summations in the surface equations of § II C 1 (although a number 
of larger fits were performed, with different initial parameter values, to make certain that 
no significant improvement was indeed possible). Finally, for a few of the "best" surfaces 
with near-optimum index ranges (and also a few very large fits), scatterplots were made of 
the error vs. the ab initio energy, and a couple of dozen contour plots of each surface were 
examined, as well as over a hundred plots of energy vs. one of the He-H 2 distances or angles 
(in these latter plots, the surfaces could be directly compared to the ab initio energies along 
some cut of the PES). Our adopted He-H 2 PES was chosen in this manner from among 
the "best" surfaces with near-optimum index ranges. A Fortran program to compute this 
analytic PES (including analytic first derivatives) is available from EPAPS 32 or from the 
authors 33 . 

We have used a relatively large number of linear parameters in order to fit our He-H 2 
data. Our final BMP surface uses 112 parameters, with a few of the non-linear parameters 
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being given reasonable constant values as discussed in § II CI, and a few others being set 
by early (small) fits and held constant thereafter (such as the parameters in and a^); in 
the final fitted surface, the 9 non-linear parameters that were refitted at that point were all 
significant (at better than the 10-cr level) , and only 2 of the 87 linear parameters had values 
smaller than their estimated errors, with only 3 others having less than a 3-cr significance 
level. These 112 parameters of our BMP surface were used to fit a total of 25 065 points 
(20 203 of which are ab initio energies), a 224:1 ratio of points to parameters (180:1 if one 
considers ab initio points only). For comparison, Truhlar and Horowitz 1 fitted 287 ab initio 
H 3 points with about 23 parameters, a 12:1 ratio. Our BKMP2 H 3 surface 4 fitted 8559 H 3 
points (7591 ab initio points) with about 120 parameters, a 71:1 ratio (63:1 for ab initio). 
Aguado et al? fitted 6101 ab initio H4 energies with 865 linear parameters, a 7:1 ratio. 
Our BMKP H4 surface 8 uses 400 parameters to fit 61547 points (48180 ab initio points), 
giving a ratio of 154:1 (120:1 for ab initio). Thus the final fitted BMP He-H 2 surface of the 
present work actually uses relatively few fitting parameters for the number of points fitted, 
minimizing the risk of "overfitting," i.e., of spurious wiggles between fitted points. 

Added confidence in our error estimates is provided by the recent extremely accurate H3 
surface of Wu et al. 5 Their ab initio energies and fitted surface are an order of magnitude 
more accurate than our earlier BKMP2 H 3 surface 4 ; they found that our H 3 error estimates 
were essentially correct, both for our H 3 ab initio energies and our fitted H 3 surface. They 
report no spurious wiggles in our BKMP2 H 3 surface with sizes larger than its quoted accu- 
racy. This suggests that our error estimates should be reasonable for the similarly-computed 
He-H2 ab initio energies of the present work, and that our fitted BMP He-H2 surface should 
likewise be free of large spurious wiggles between fitted points. Finally, random "interstitial" 
unfitted conformations were also tested directly, as described in § III A 1. 

III. DISCUSSION 
A. Accuracy of analytic He-H2 surfaces 

Our adopted surface (the "BMP" He-H 2 PES) contains a total of 112 parameters; even 
much larger surfaces (up to 175 parameters) yielded no significant improvement (having rms 
values only a few percent smaller). Figure 5 shows a scatterplot of the deviations relative to 
our ab initio energies. The surface fits best at low energy (E < 0.15 E h) i.e., energies below 
about twice the H 2 -molecule dissociation energy, relative to the energy of an equilibrium 
H 2 molecule plus a distant He atom). However, even up to quite high energies, the surface 
lies within a few mE h of most of the ab initio points; it turns out that most of the outliers 
are due to the fact that the surface can only produce an approximate fit to the conical 
intersection with the first excited state (this is discussed in more detail below, in § III A 3). 

Table II compares the rms errors of four earlier analytic He-H 2 surfaces with rms errors 
of the adopted (BMP) surface of this work, for various subsets of ab initio and generated 
energies; the modified Muchnick & Russek (modMR) surface of § II D 2 (which was used to 
generate some of the constraining van der Waals points for non-equilibrium H 2 ) is also shown. 
"Energy- weighted" deviations were used to obtain the rms values of Table II, i.e., energies 
above about two H 2 -dissociation energies were given reduced weight we according to the 
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formula of Equation (22). The only exception to this is the last line of Table II (subset 44), 
where the fully weighted rms values are reported for all points used in our He-H 2 fit (for our 
BMP surface, this is the rms value that the fitting program tried to minimize, as described 
in § IID1 and § IID3). 

It should be noted that the earlier analytic He-H 2 surfaces were fitted to ab initio points 
in relatively restricted regions, compared to the region covered by the ab initio and generated 
points of the present work. The rms values given in italics in Table II indicate subsets of 
points that cover regions where the corresponding surface was not fitted, i.e., where energies 
predicted by these previous surfaces are being tested by our new ab initio energies. (The 
only unconstrained region of our fitted BMP surface is at very small interatomic separations, 
i.e., very high energies; this is discussed in § III A 6.) 

The Wilson, Kapral, & Burns 17 (WKB) surface was fitted to a set of ab initio energies 
in the region 0.75 a < r < 5.0 a , 0.0 a < R < 5.0 a , plus the three pair-wise potentials. 
The Dove & Raynor 20 (DR) surface was essentially the same, with improved versions of the 
pair-wise potentials. Since the interaction part of these surfaces was fitted using only four 
free parameters, they lack sufficient flexibility for accuracy even in their fitting region. 

The Schaefer & Kohler 22 (SK) surface was fitted to a set of ab initio energies in the 
region 1.28 a < r < 1.618 a , 3.0 a < R < 15.0 a plus a somewhat less accurate set in 
the region 0.9 a Q < r < 2.0 a , 1.5 a < R < 8.0 a ; it is quite accurate in these regions. 
However, it is designed to be interpolated on a grid in the region 0.9 a < r < 2.0 a , 
R > 1.6 a , and therefore can be highly inaccurate when extrapolated outside this region 
(just how inaccurate depends on the extrapolation formula used — we just used an spline 
in r, even for extrapolation, but fitted exponentials to extrapolate the R spline to small R). 

The Muchnick & Russek 24 (MR) surface was fitted to a set of energies in the region 
1.2 a < r < 1.6 a , 0.0 a < R < 15.0 a ; the very similar modMR surface of § IID2 
differs from the MR surface only in the parameters describing the H2 — He and H — He van 
der Waals wells. These surfaces are accurate in the part of their fitting regions that lies at 
R ^ 3 a . However, although the MR surface was fitted to ab initio energies down to R = 0, 
it is a relatively poor fit at small R (subset 9 in Table II). This is due both to the fact that 
Muchnick & Russek 24 had only a few ab initio points in this region and to the fact that a 
good fit in this region requires a good deal of flexibility in the fitted surface (we found in 
our first, smallest fits that significant improvements at small R could be obtained only at 
the expense of worsening the fit in the van der Waals well, unless more parameters were 
added). Even at intermediate R (2 a < R < 3 a : subsets 6, 10, and 14 in Table II), the 
SK surface does somewhat better than the MR surface, provided that r is not too large, i.e., 
r < 2 a (our BMP surface is of course significantly better than either, in these regions of 
small-to-intermediate R) . 

Table II demonstrates that, overall, the MR surface 24 is the best of the previous surfaces 
considered here. Nonetheless, even the MR surface does quite poorly overall, with an rms 
error of 28.6 mE h relative to our ab initio energies (subset 1 in Table II). The present BMP 
surface is an order of magnitude improvement, with an rms error of 1.42 mE h relative to the 
ab initio energies — of the same order as the estimated ab initio error of ~ 0.6 mEh- The 
advantage of the BMP surface is largest in the region that can be sampled by dissociative 
collisions, but even below the dissociation energy (i.e., E < 0.0 E^. subset 3 in Table II) 
the present BMP surface (rms of 0.48 mE h ) remains a significant improvement over the MR 



19 



surface (rms of 4.48 mEh). 

For the restricted region of not-too-large H 2 -molecules (r < 2 a ) plus a not-too-close 
He atom (R > 3 a ), subsets 7, 11, 15, and 31 in Table II show that both the SK and MR 
surfaces 22 ' 24 do nearly as well as the present BMP surface, but the earlier WKB and DR 
surfaces 17 ' 20 are an order of magnitude worse even in this region (and have no van der Waals 
well at all). As one might expect, both the SK and MR surfaces do better when extrapolated 
to small r (subsets 4 to 7 in Table II) than when extrapolated to large r (subsets 16 to 19), 
due presumably to the fact that reducing r tends to reduce the anisotropy with respect to 7 
at reasonable R values, while for large r the anisotropy can become very large (in fact, for 
R ~ r/2, the He atom approaches one of the H atoms for 7 — > 0). 

We performed a number of further tests to assess in more detail the quality of our BMP 
surface (and of the previous surfaces), as discussed in the following subsections. 

1. Test for "interstitial" wiggles in the BMP surface 

The rms errors of the 3500 unfitted "random" conformations of § II A 3 (subsets 23 to 25 
in Table II) were compared to the rms errors of the 3500 fitted "random" conformations of 
§ II A 2 (subsets 20 to 22), which had been generated with the same probability distribution. 
The similarity of these two sets of rms values is evidence that our BMP surface does not 
have large spurious wiggles between fitted ab initio points. The fact that the rms values 
for these sets of "random" conformations are also similar to the overall ab initio rms values 
(subsets 1 to 3, dominated by the main grid of § II A 1) yields added reassurance. 

These "random" points of these subsets (20 to 25 in Table II) avoid duplicate conforma- 
tions by a relatively small "avoidance margin" of ~ 0.01 a (see § II A 2), in order to keep 
the fraction of "avoided" conformation space small even in regions where the main grid of 
conformations was dense. We also considered rms errors of our BMP surface for "maximally 
interstitial" cases with larger avoidance margins: e.g., for a case where points from fitted 
subset 20 were discarded if they lay within a 3% avoidance margin, the remaining 2949 
points had an (energy- weighted) rms error of 1.356 mEh (0.528 mEh for the 1140 points 
with E < 0), while the corresponding 2742 points from unfitted subset 23 had an rms error 
of 1.374 mEh (0.536 mEh for the 1063 points with E < 0). (Recall that the typical main-grid 
spacing was ~ 7% in r and R A , ~ 7% in 9 for 9 > 120°, and ~ 13% in 9 for 9 < 120°: see 
§ II A 1.) As expected, there was good agreement between rms values of fitted and unfitted 
"random" subsets independent of the size of the avoidance margin. 

It would of course have been very surprising if any interstitial wiggles had showed up; 
as discussed at the end of § II D 4, the density of fitted ab initio points is high enough that 
there is no room for interstitial wiggles (unless we had used a great many more parameters 
in the fit). However, there are nonetheless certain regions where the BMP surface does have 
errors significantly larger than its typical ones, as discussed further below. 

2. Contour plots of the BMP surface 

The contour plots of Figure 6 show the overall shape of our fitted BMP surface for four 
orientations 7 of the H 2 — He separation R relative to the H — H separation r. Note that 
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the total potential VHeH 2 is represented in Figures 6 and 7 (not the interaction energy Vi n t), 
and thus its gradient gives the total force on the H and He atoms. 

In Figure 6a, with 7 = 0°, the He atom lies very close to one of the H atoms for R pa 0.5 r, 
yielding a double-sided "wall" in the contour plot. In Figure 6b, this "wall" is replaced by a 
ridge, since for 7 = 30° the He atom is closest to the H atom when R = 0.5 r. The ridge-like 
structure at lower left in Figure 6c (7 = 60°) is probably due to the attempt by the surface 
to fit the nearby conical intersection with the first excited state, which is discussed in more 
detail in § III A 3 below. 

From the slopes of the contour lines in Figure 6d (a T-shaped orientation, i.e., 7 = 90°), 
one can see that, as R decreases, the He atom does not begin to "push apart" the H atoms 
until it is quite close to the H2 molecule: the partial derivative {dV/dr)n n remains positive 
for r > 1.5 a until R reaches ~ 1.5 a , between the 100 and 200 mEh contours. (There is 
a large repulsive force between the He atom and the H 2 molecule that begins much further 
out, as may be seen from the significant reduction in depth of the H 2 potential well as R is 
reduced below ~ 4 o ). The cuts at constant r in Figure 7 also illustrate this effect — the 
curves for larger r values lie above those for r = 1.4 a , except at small R. This illustrates an 
effective barrier to dissociation of H 2 by He: at least in the T-shaped geometry, the He atom 
must have a high enough energy to sample the part of the surface where (dV/dr) < 0, i.e., 
considerably more than just the the H 2 dissociation energy. 

For an equilibrium H 2 (r = 1.4 a ) and 7 = 90°, Muchnick & Russek 24 found that, 
with their surface, the presence of the He atom actually caused an inward force on the 
H 2 molecule for 2.36 a < R ^ 6.5 a , i.e., (dVi nt /dr)R n > in this region; for the BMP 
surface of the present work, the corresponding region is 2.03 a < R < 6.29 a . As Muchnick 
& Russek 24 pointed out, this "force reversal" is a real effect, even visible in the ab initio 
energies. Nonetheless, this "force reversal" is only in the interaction energy, and is a minor 
effect compared to the stability provided by H 2 -molecule potential; as may be seen from 
Figure 6, the equilibrium H 2 size is almost independent of R, except at small R. 

Figure 7 also illustrates that the nearly-vertical parts of the contours at the lower left 
in Figure 6d (7 = 90°) are a real reflection of the ab initio data (and not an artifact in the 
fit arising from fitting the smaller-7 region where the He atom passes close to one of the 
H atoms). For 7 = 90° and r ~ 3 a , there really is a (well-fit) plateau in the energy as a 
function of R for R < 1.5 a . For the r < 2.55 plotted in Figure 7, this plateau is 

more tilted and lies at higher energy; the "wiggles" of up to ~ 10 mEh visible in the fitted 
BMP surface there are probably due partly to the reduced weight given to high-energy 
points, but are probably also connected to the relatively nearby conical intersection with 
the first excited state (see § III A 3 below). 

Note that in Figure 6 the ~ 0.3 mEh "basin" near {r = 10.5 a , R = 0}, and also the 
~ 0.3 mEh "ridge" in Figures 6c and 6d near {r = 12 a„, R = 2 a,,}, are probably spurious 
features, where our fitted BMP surface overcompensates slightly due to having to fit the 
nearby repulsive walls (the spurious features being comparable in size to the expected errors 
in the ab initio energies in this region). Figure 8 illustrates that there is no real evidence for 
this "basin" and "ridge" in the ab initio energies or generated points. These (presumably) 
spurious features in our BMP fit are smaller even than the rms error with which the surface 
fits the low-energy points, and furthermore lie in a region of the surface not likely to be 
sampled, i.e., by the time the H atoms are that far apart in a dissociative collision, the 



21 



He atom is not likely to lie half-way between them. Consequently, there seemed to be little 
point in attempting to smooth out the "ridge" and "basin" by increasing the weight on the 
points in this region (especially since this might lead to a poorer fit elsewhere). 

Figure 8 also shows that the previous He-H 2 surfaces have a larger spurious "bulge" or 
"basin" than our BMP surface in this region. Through the use of A d (r) from Equation (13), 
which has an exponential cutoff at large r, the modMR surface of § II D 2 largely eliminates 
the spurious "basin," but, like the MR surface itself, is a poor fit to the base of the repulsive 
"wall" in the 7 = 0° direction. (In fact, the limit r — > oo at constant Ra was not constrained 
at all in the Muchnick & Russek 24 fit, and thus it is not surprising that their resulting MR 
surface and the similar modMR surface yield a very poor fit to the pair- wise H-He potential.) 
The WKB and DR surfaces did include a pair-wise H-He potential, yielding a reasonable 
fit to the base of the repulsive H-He "wall" in the 7 = 0° direction (although their H-He 
potentials were less accurate at higher energies). However, their three-body terms did not 
have a van der Waals well, but rather a relatively long positive "tail" at large separations, 
resulting in the "bulge" (in the 7 = 90° direction) shown by these potentials at small R in 
Figure 8 (where r = 10.6 a and Ra = Rb > 5.3 a ). 

In the following subsections, we consider the quality of our BMP fit in various other 
regions of He-H 2 conformation space. 

3. Conical intersection with the first excited state 

The conical intersection of the ground state with the first (electronic) excited state forms 
a curved (1-dimensional) line in the 3-dimensional conformation space of He-H 2 . Since we 
used the MRD-CI program to calculate the first few excited states as well as the ground 
state, we were able to estimate the position of this conical intersection by considering the 
energy difference between the ground state and first excited state in the ab initio energies. 
Figure 9 shows this estimated position: the squares, with ground state and first excited 
state lying at nearly the same energy, should lie closest to the conical intersection. The 
position of the conical intersection is thus determined to an accuracy of ~ 0.03 a in Ra 
and ~ 3° in 0, for 1 a < r < 6 a . Figure 9d shows that, in contrast to the H 3 surface, the 
conical intersection for He-H 2 does not in general occur at high-symmetry conformations 
(i.e., conformations having 7 = 90°). 

Somewhat better estimates of the position of the conical intersection (as well as estimates 
of its opening angle in the coordinate directions orthogonal to the line of the conical inter- 
section, as a function of position along that line) could be obtained by plotting the energies 
of the ground and lowest excited state(s) as a function of position for gridpoints near the 
estimated conical intersection. Even better accuracy could be obtained by the computation 
of new ab initio energies to home in on the position of the conical intersection (e.g., as 
was done for H 4 in a limited region by Boothroyd et al. 8 ). Procedures such as these would 
be required if one wished to include the conical intersection explicitly in a fitted surface; 
however, they are beyond the scope of the present work. 

As may be seen from Figure 9e, the conical intersection occurs at energies E > 0.2 E h: 
i.e., above twice the H 2 dissociation energy relative to H 2 + He. Figure 10 shows various 
cuts through points on the conical intersection (note that Figs. 10a, 10b, and 10c show 
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three different cuts through the same conical intersection point, while Fig. lOd shows a 
fourth cut through two different conical intersection points). The BMP surface can be seen 
to "round off" the conical intersection, missing by as much as 20 mE h (although it still 
does much better than the previous surfaces in this region of the He-H 2 surface). Figure 9f, 
showing errors of the BMP surface relative to the ab initio energies at points near the conical 
intersection, also shows the surface lying quite far below the cusp of the conical intersection 
(black squares) for r > 1 a , and (at least in some cases) overcompensating at points slightly 
further away (black dots). As a result, the conical intersection is the source of many, perhaps 
most, of the outliers at E > 0.2 E h in the scatterplot of Figure 5. 

4- The Hi + He van der Waals well 

Figure 11 compares fitted He-H 2 surfaces to ab initio van der Waals energies — note 
that the unusual "shifted-logarithmic" energy scale emphasizes effects at the bottom of 
the van der Waals well, while also showing a significant portion of the "repulsive wall" at 
smaller R. One can see that our BMP He-H2 surface fits the depth of the van der Waals well 
to within about 1 fj,Eh, when compared to the accurate ab initio energies of Tao 31 (which 
have 1.28 a < r < 1.618 a ); from subset 31 in Table II, the rms difference is 0.8 /mE^ over 
the range 6 a < R < 10 a . The previous SK and MR surfaces 22 ' 24 agree nearly as well 
(i.e., within a few fiE^, rms values of 1.8 and 2.6 fiE^, respectively), since they were fitted 
to the fairly accurate ab initio energies of Meyer et al. 29 However, the earlier WKB and DR 
surfaces 17 ' 20 do not contain a van der Waals well at all — they are visible only at the upper 
right of Figure 11. 

Note that the "torque reversal" at R ~ 6.1 a noted by Muchnick & Russek 24 is also 
visible in Figure 11. For small angles 7, the van der Waals well is slightly deeper but 
the "repulsive wall" is slightly further out than for 7 near 90°. Thus the 7 = 0° (linear) 
orientation is energetically favored at R > 6.1 a , while the 7 = 90° (T-shaped) orientation 
is favored at smaller distances R. 

Figure 12 considers the variation of the van der Waals well with the size r of the H 2 
molecule (using the same energy scale as Fig. 11). For r = 1.28, 1.449, and 1.618 a , 
accurate ab initio energies from Tao 31 are available, which our BMP surface fits very well 
(rms of 0.8 fiE h : see subset 31 in Table II). For the more extreme distances r = 0.9, 1.1, 
1.8, and 2.0 a , only the much-less-accurate "gen-SK vdw" and "gen-MR vdW" generated 
points of § II A 4 were available. The former were obtained from the SK surface 22 , which at 
those r values was fitted to the lower-accuracy extreme-position ab initio energies of Meyer 
et al. 29 ; the latter were obtained from the modMR surface of § II D 2, fitted to the accurate 
near-H 2 -equilibrium ab initio energies from Tao 31 with the r-dependence taken from Thakkar 
et al. 23 These two sets of generated points agree with each other to within about 5 [iE h in 
the van der Waals well, and are fitted to roughly this accuracy by our BMP surface (rms 
of ~ 3 fiEh over the range 6 a < R < 10 a , from subsets 30 and 32 in Table II; largest 
difference ~ 10 /iEh). As the size r of the H 2 molecule is increased, the van der Waals well 
moves slightly further out; the depth also changes slightly, but this is less well constrained 
by the ab initio energies. 

At r > 2 a , there are no accurate ab initio energies available in the region of the van der 
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Waals well (some of our ab initio energies may lie nearby, but these have estimated errors 
larger than the typical depth of the van der Waals well). The generated "vdW" energies 
in this region are thus expected to be quite inaccurate — the SK and modMR surfaces 
differ from each other by ~ 0.1 mE h in this region. These generated points were therefore 
given quite low weight, only enough to keep the fitted surface from developing large spurious 
features in this region. This is the cause of the relatively large rms error of 0.1233 mEh of 
the BMP surface shown by Table II (subset 33) for the "vdW" at r = 2.4,3.0,4.0 a and 
6 a < R < 10 a . 

5. The interaction region 

As mentioned above, Table II shows that the earlier SK and MR surfaces 22 ' 24 do relatively 
well for not-too-large H 2 -molecules (r < 2 a ), even up to relatively large interaction energies. 
An example of this, with r = 1.4 a , is illustrated in Figure 13a (note the logarithmic energy 
scale). Even at R ~ 1.3 to 2 a , these two surfaces have errors of only a few percent for 
near-equilibrium H2 — the BMP surface of the present work is an order of magnitude more 
accurate, with errors of a fraction of a percent in this region, but nonetheless all three of 
these surfaces are almost indistinguishable in Figure 13a. Even the much-earlier WKB and 
DR surfaces 17,20 are not so very much worse in the outer part of this region, but, as may be 
seen in Figure 13a, they have an unphysical "hole" at R ^ 1.5 a for small angles 7 (i.e., 
the "hole" is at small H — He separations Ra)- 

Table II shows that, for all the fitted surfaces considered there, errors at small r (i.e., 
r < 1.2 a : subsets 4 to 7) are very similar to the errors at r ~ 1.4 a (subsets 8 to 11) that 
are illustrated in Figure 13a. In contrast, at larger H — H separations r (such as would arise 
from dissociative collisions and highly-excited H 2 molecules), the SK surface 22 is undefined, 
and even the MR surface 24 does quite poorly, with errors of order 50% visible in Figure 13b, 
c, and d (for r = 4 a ). At these larger r values, the WKB and DR surfaces 17,20 have errors 
of order a factor of 2, and Figure 13b shows that the "hole" in these surfaces appears at 
lower interaction energy (and, in fact, at a lower total energy, though this is less obvious 
from the figure). A systematic numerical search found that the lowest part of the lip of this 
"hole" lies at energies E = 53 and 87 mEh for the WKB and DR surfaces, respectively - 
less than two H 2 dissociation energies above equilibrium H 2 + He. In this region, the BMP 
surface of the present work is much more accurate than any previous surface, with typical 
errors only slightly larger than for equilibrium H 2 molecules. This is also illustrated by the 
other cuts through this region shown in the previous Figure 10. 

6. Extrapolation to very short distances 

The SK surface 22 contains no repulsive Coulomb terms, and thus does not have the 
correct form at small H — He separations — an example can be seen in Figure 13a. The 
WKB and DR surfaces 17,20 do even worse: they turn over and become negative at small 
H — He separations (Ra or Rb 1 a : the "hole" in these surfaces), as may be seen in 
Figure 13a and b. The MR surface 24 behaves roughly as an exponentially damped Coulomb 
repulsion for small Ra, yielding reasonable behavior at short distances. 
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No constraints were placed on our fitted BMP surface at H — He distances less than 
those of the most compact ab initio geometries (namely Ra = 0.6 a ), which lie at energies 
E > 1 E h . Some of our earlier fitted surfaces turned over into an unphysical "hole" at 
shorter distances than this, but our final BMP He-H 2 surface becomes strongly repulsive 
there instead. An example of this behavior may be seen in Figure 13a and b: above ~ 1 Eh 
the BMP curve (solid line) rises very steeply. Extrapolating from the ab initio energies, 
the short-distance behavior at H — He distances Ra < 0.6 a ought to be intermediate 
between the BMP curve and that of the MR surface 24 (long-dashed line), but closer to 
the latter for R A ^ 0.55 a . For cases where an extremely hard repulsive core at Ra < 
0.6 a is undesirable, the Fortran program of our BMP analytic surface 32 ' 33 contains an 
option allowing the user to switch over at small Ra or Rb (gradually, with continuous first 
derivatives) to alternate forms, which behave more like an exponentially damped Coulomb 
repulsion there. 

All of the analytic surfaces include the H 2 -molecule potential V H2 as a separate term in the 
surface formulae. The WKB and MR papers 17 ' 24 left the choice of V H2 formula unspecified, 
while the Fortran program for the DR surface 20 used the spline fit of Kolos & Wolniewicz 48 
for Vh 2 (which is almost as good as Schwenke's 28 Vh 2 formula). For consistency, we used 
the accurate H 2 -molecule potential of Schwenke 28 for all the surfaces. This formula behaves 
well at most H — H separations r, but begins to behave poorly at r < 0.1 a . Our Fortran 
program therefore switches over at r ?a 0.12 a to a polynomial in 1/r, which in turn becomes 
a pure 1/r repulsive potential for r < 0.045 a . (The first derivative with respect to r is 
continuous at these switchover points r ~ 0.12 a and r 0.045 a .) 

B. Prospects for further improvement 

Significant improvements over the BMP surface of this paper are possible, but would 
require a major effort. 

Improvements in the van der Waals well for strongly non-equilibrium H 2 molecules (r < 
1.2 a and r > 1.6 a ) would require very accurate ab initio calculations in this region (such 
as those performed for equilibrium H 2 molecule sizes by Tao 31 ). Once these ab initio energies 
were available, increased flexibility (i.e., larger index ranges) in the Va and Vd terms should 
allow the van der Waals well to be well represented over a wide range of H 2 molecule sizes. 

Improvements in the interaction region would require, in the analytic functional form, 
explicit inclusion of the conical intersection with the first excited state. To do this, one would 
first have to pin down its location and its opening angle in the perpendicular directions (it is 
possible that these could be determined sufficiently well from the existing ab initio energies 
of the present work; the alternative would be to use an explicit search via new ab initio 
computations, as was done for H4 in a limited region by Boothroyd et al. 8 ). The position 
and opening angle would then have to be fitted to an analytic form, which could then be 
added to the fitted surface — the magnitude of this new term (and perhaps its extent 
in the directions perpendicular to the locus of the conical intersection) would need to be 
parameterized and fitted as a function of position along the locus of the conical intersection. 
The other terms in the fitted surface would then not have to try to fit the sharp peak very 
near the conical intersection, allowing a considerable improvement in accuracy there. 
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IV. CONCLUSIONS 



A new set of 23 703 ab initio energies was computed for He-H 2 geometries where the 
interaction energy is expected to be non-negligible, using Buenker's multiple reference (single 
and) double excitation configuration interaction (MRD-CI) program 40 ' 39 ; the lowest excited 
states were computed as well as the ground state energy. These new ab initio energies 
have an estimated rms "random" error of ~ 0.2 mEh and a systematic error of ~ 0.6 mEh 
(0.4 kcal/mol). The position (in the 3-dimensional conformation space of He-H 2 ) of the 
conical intersection between the ground state and the first excited state has been roughly 
mapped out; unlike H 3 , this conical intersection for He-H 2 does not lie at high-symmetry 
conformations, but rather along a curved line in conformation space. 

These new ab initio energies were used to test previous analytic He-H 2 surfaces. Even 
the best of the previous surfaces, that of Muchnick & Russek 24 , does quite poorly for very 
large H 2 -molecule sizes (such as would be encountered in dissociative collisions or highly 
excited H 2 molecules), although both this surface and that of Schaefer & Kohler 22 are fairly 
accurate for not-too-large H 2 molecules (r < 2 a ) with a not-too-close He atom (R > 3 a ). 

A new analytic He-H 2 surface (the BMP surface) was fitted to 20 203 of the new ab initio 
energies (and to an additional 4862 points generated at large separations) - - the other 
3500 new ab initio energies were used only to test "interstitial" conformations. This BMP 
surface fits the interaction region at the "chemical accuracy" level (~ 1 mE h ) required for 
reaction dynamics; the overall energy-weighted rms error of 1.42 mEh (0.89 kcal/mole) is 
comparable to the accuracy of the ab initio energies (note that this energy- weighting specifies 
lower weight for high-energy points, i.e., E > 0.2 Eh, but a weight of unity everywhere else: 
see § II D 1). For the 14 585 ab initio energies that lie below twice the H 2 dissociation energy, 
the new BMP He-H 2 surface has an rms error of 0.95 mE h (0.60 kcal/mole). This surface 
is an order of magnitude better than previous surfaces in the interaction region, and also 
yields a slight improvement in the fit to the recent van der Waals energies of Tao 31 . 

For relatively compact conformations (i.e., with energies above twice the H 2 dissociation 
energy), the conical intersection between the ground state and the first excited state is the 
largest source of error in the analytic surface. The BMP surface "rounds off" the conical 
intersection, yielding errors of up to ~ 20 mEh there. The position of this conical intersection 
forms a curved 1-dimensional locus in the 3-dimensional conformation space of He-H 2 ; its 
approximate position has been mapped out, but trying to include the conical intersection 
explicitly in an analytic surface would require a more accurate map of its position, and is 
beyond the scope of the present paper. 

The ab initio energies and a Fortran program for the analytic BMP He-H 2 surface of the 
present work are available from EPAPS 32 or from the authors 33 . 
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TABLES 





TABLE I. 


Gaussian basis sets 






Tvnp 




(HopflR PI PTlt 

\J \J Cll 1 \j 1C111 


H: (9s3p2d)/[4:s3p2d] 


S 


887.22 


0.000112 






123.524 


0.000895 






27.7042 


0.004737 






7.82599 


0.019518 






2.56504 


0.065862 






0.938258 


0.178008 




s 


0.372145 


1.0 




s 


0.155838 


1.0 




s 


0.066180 


1.0 




p 


2 11 75 


1 




V 


77 


i n 




p 




1 n 




Uj 


1 76 


i n 




A 

Hi 


62 


1 


He: (Ws4p2d)/[5s4p2d] 


S 


3293.694 


0.00010 






488.8941 


0.00076 






108.7723 


0.00412 






30.17990 


0.01721 






9.789053 


0.05709 






3.522261 


0.14909 




s 


1.352436 


1.0 




s 


0.552610 


1.0 




s 


0.240920 


1.0 




s 


0.107951 


1.0 




p 


5.823125 


1.0 




p 


2.1175 


1.0 




p 


0.77 


1.0 




p 


0.28 


1.0 




d 


2.83871 


1.0 




d 


1.0 


1.0 
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TABLE II. Energy-weighted a rms errors of analytic He-H2 surfaces, in mEh 



Surface I.D.: b 




WKrJ 


I ~\ T 1 

DK 


SK 


MR 


modMR 


I -) TV T 1 ~l 

BMP 


No. of parameters: 




4 


4 


[750J 


19 


27 


112 


Subset 




rms 


rms 


rms 


rms 


rms 


rms 


1 




All fitted ab initio 


20 203 


100.90 


100.90 


5.1E+5 


28.60 


28.35 


1.417 


2 




... E < 0.174 E h A 


14 585 


31.80 


40.40 


6.0E+5 


20.24 


19.89 


0.955 


3 




... E < 0.0 E h e 


7177 


11.62 


10.58 


1.IE+4 


448 


444 


0.481 


4 




... 0.5 a < r < 1.2 a 


3744 


26.75 


20.95 


21.35 


1.80 


1.80 


0.875 


5 




R<2a 


1646 


37.73 


29.57 


32.18 


10.81 


10.81 


1.205 


6 




2 a < R < 3 a 


857 


17.89 


14.18 


1.323 


6.20 


6.21 


0.624 


7 




R > 3 a 


1241 


7.11 


4.97 


0.449 


0.418 


0.421 


0.336 


8 




... 1.2 a <r< 1.62 a 


2311 


48.62 


47.82 


22.14 


15.79 


15.79 


1.003 


9 




It ^ Z Bo 


907 


75.84 


74.80 


35.31 


24.69 


24.68 


1.463 


1(J 


2 a < i? < 3 a 


575 


18.72 


18.04 


1.497 


6.39 


6.40 


0.690 


1 1 
11 


R > 3 a 


829 


7.39 


5.19 


0.433 


0.445 


0.450 


0.365 


iz 


... 1.62 a < r < 2 a 


1918 


74.01 


76.16 


116.11 


25.34 


25.33 


1.351 


13 


R<2a a 


735 


117.85 


121.07 


181.64 


39.59 


39.58 


2.068 


14 


2 a < R < 3 a 


518 


22.31 


25.29 


3.42 


12.36 


12.34 


0.723 




i? > 3 a 


665 


7.71 


5.51 


0.458 


0.865 


0.889 


0.358 


16 


... 2 a < r < 4 a 


T1 OCX 

I lob 


130.2 f 


133.13 


i i TP t 1 

l.l£j+4 


3 1.63 


3 1. 64 


1. 1 9 ( 


17 


it ^ O U/Q 


2806 


10.16 


8.30 


24.15 


1.64 


1 11 


0.456 


18 




5094 




109. 56 


1.0E+6 


29.32 


28.34 


1.318 


19 




1244 


8.44 


6.40 


84.99 


18. 41 


11.81 


0.588 


20 


Fitted "random" ab initio 1 


3500 


51.82 


64.25 


2.9E+5 


21.36 


26. 95 


1.411 


21 


... E < 0.174 E h d 


2962 


36. 35 


39. 61 


3.1E+5 


20. 81 


20.32 


0.994 


22 


... E < 0.0 E h e 


1367 


13. 92 


13. 03 


5.11 


5.01 


5.04 


0.527 


23 


Unfitted "random" ab initio^ 


3500 


51.01 


63.35 


2.6E+5 


21.14 


21.31 


1.400 


24 


... E < 0.174 E h d 


2943 


35.58 


38.41 


2.8E+5 


21.05 


20.52 


0.964 


25 


... E < 0.0 E h e 


1349 


14.10 


13.12 


14.50 


4.93 


4.89 


0.538 


26 


All "vdW" 


2145 


2.031 


1.698 


0.134 


0.609 


0.609 


0.129 


27 


... R < 6 a 


587 


3.61 


3.03 


0.251 


1.164 


1.164 


0.237 


28 


... 6 a < R < 10 a 


1098 


0.930 


0.830 


0.0319 


[0.0114] 


[0.0013] 


0.0507 


29 


r = 0.6,0.7,0.8 a h 


186 


0.491 


0.404 


0.0023 


[0.0020] 


[0.0000] 


0.0031 


30 


r = 0.9,1.1 ao 1 


248 


0.131 


0.640 


[0.0011] 


[0.0014] 


[0.0011] 


0.0017 


31 


r = 1.28,1.449,1.618 a Q j 


231 


1.004 


0.906 


0.0018 


0.0026 


0.0011 


0.0008 


32 


r = 1.8,2.0 ao 1 


248 


1.139 


1.038 


[0.0024] 


[0.0042] 


[0.0024] 


0.0037 


33 


r = 2.4,3.0,4.0 a h 


185 


1.084 


0.952 


0.0922 


[0.0271] 


[0.0000] 


0.1233 


34 


... i? > 10 a 


460 


0.213 


0.210 


0.0004 


0.0005 


0.0001 


0.0001 


35 


All "H-He" 


1887 


114.69 


106.25 


4.9E+4 


44.54 


44.00 


1.122 


36 


... i?A < 2 a 


1008 


156.90 


145.35 


1.2E+3 


59.88 


59.18 


1.455 


37 


... 2 a < Ra < 4 a 


468 


3.65 


2.75 


4.5E+3 


16.59 


16.19 


0.714 



30 



38: ... 4 a < Ra < 6 a 
39- i?4 > 6 a n 


228 
183 


0.416 
0.0431 


0.0959 
0008 


3.8E+3 
1.6E+5 


0. 680 
0068 


0. 639 
0038 


0.0978 
0.0014 


A C\ A 11 ££TT TT i TT*1 

40: All H-He + H 


830 


26. 03 


22.19 


7.6E+6 


19. 82 


19.34 


0.874 


41: ... R A < 3 a 


141 


63.15 


53.84 


2.4E+4 


48.03 


46.88 


2.096 


42: ... 3 a < i?^ < 6 a 


201 


0.950 


0.602 


8.9E+5 


2.02 


1.63 


0.266 


43: ... i? A > 6 a 


488 


0.184 


0.156 


9.9E+6 


0.0813 


0.0051 


0.0341 


44: All fitted points; full-weight a 


25 065 


329.02 


310.15 


5.2E+7 


28.70 


28.41 


1.383 



a For energy- weighted rms errors, deviations get only weight we from equation (22) (unity for 
E < 0.2 Eh, reduced weight at higher energy), as opposed to full-weight rms errors in last line 
of table (subset 44). Note that the rms value for a surface is italicized for subsets containing 
conformations in regions where that surface was not fitted. 

b WKB: Wilson, Kapral, & Burns (Ref. 17), DR: Dove & Raynor (Ref. 20), SK: Schaefer & Kohler 
(Ref. 22), MR: Muchnick & Russek (Ref. 24), modMR: slightly modified Muchnick & Russek 
surface of § II D 2, BMP: adopted surface of the present work. 

c Schaefer & Kohler present 3 Legendre coefficients at 48 i?-values for each of 5 r-values (namely, 
r = 0.9, 1.28, 1.449, 1.618, and 2.0 a ), plus formulae for very large R. 

d E < 0.174 Eh corresponds to points lying below about twice the H2 dissociation energy, relative 
to a H2 molecule plus a distant He atom. 

c E < 0.0 Eh corresponds to points lying below the H2 dissociation energy, relative to a H2 molecule 
plus a distant He atom. 

f These subsets 20 to 22 are also contained in subsets 1 to 3, respectively. 

g These subsets 23 to 25 were obtained using the same probability distribution as subsets 20 to 22, 
but were not included in the actual fitting process for the BMP surface of the present work. 
h The accuracy of these generated points is unknown. They were intended only to yield "reasonable" 
values and prevent wild excursions in the fitted surfaces, and were given relatively low weight in 
the fit, resulting in the relatively large "BMP" rms for these points at large r. The "modMR" 
surface was used to generate them, and thus by definition has zero rms here (likewise, the very 
similar "MR" surface has a small rms here) ; this is indicated by the square brackets enclosing their 
rms values. 

1 The "SK" and "modMR" surfaces were used to generate these points, and thus they (and the 
"MR" surface) have small rms values here (again indicated by square brackets), 
j These points comprise the ab initio energies of Tao (Ref. 31) and points generated from them, 
and thus these rms values provide a measure of the accuracy for all the surfaces in the van der 
Waals well. 
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H atom position: for r = 2.35 a 
He atom: ab initio main grid, at r = 2.35 a 
He atom: random ab initio, 2.275 a < r < 2.45 a 
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H atom position: for r = 2.35 a 
He atom: ab initio main grid, at r = 2.35 a 
He atom: van der Waals H 2 + He, at r = 2.4 a 



HS..S..S..S. 



1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 

position y 3 of atom (a ) 



position y 3 of atom (a ) 

FIG. 2. The grid of fitted points for a given H2-molecule size r, showing positions of the 
He atom relative to a molecule lying on the z-axis and centered at the origin, (a) The main grid 
of ab initio points (solid dots), for r = 2.35 a (the H-atom position is shown by the open circle); 
crosses show positions of "nearby" random ab initio points (in the same "r-bin", namely, with 
2.275 a < r < 2.45 a ). (b) The grid for the van der Waals H2 + He points (open squares), for 
r = 2.4 a ; the main grid of part (a) is also visible (solid dots). 



33 



29 
28 
27 
26 
25 
24 
23 
22 
21 

_ 20 
o 

ns. 19 
E 18 

O 17 



CO 



N . 



° 13 

OJ 

N 12 

.2 11 

w 10 
o 

Q- 9 

8 
7 
6 
5 
4 
3 
2 
1 




H atom position: for r = 9.25 a 
He atom: ab initio main grid, at r = 9.25 a 
He atom: "approximate H-He + H", at r = 9.25 a 
He atom: "accurate H-He", 9.025 a < r < 9.475 a 
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position y 3 of atom (a ) 



FIG. 3. The grid of fitted points for a large H2 separation r = 9.25 a ; unlike Fig. 2, the ab initio 
points (solid dots) all lie relatively near the origin, between the H atoms. The "approximate 
H-He + H" grid (open squares) at r = 9.25 a was designed to avoid other grids, including any 
"nearby" points on the "accurate H-He" grid (crosses — with 9.025 a < r < 9.475 a ). 
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FIG. 4. Our fits Ad(r) and a^(r) (solid and dashed lines, respectively) to the asymptotic pa- 
rameters Ad and ad of Thakkar et al. (Ref. 23: circles and squares), as a function of the H2 molecule 
size r (note that actual ad values have been multiplied by a factor of 30 so as to be visible on the 
same plot as Ad). Values extrapolated to r < 0.8 a are of little practical relevance, as H2 is un- 
bound at the corresponding energies. The linear forms of Aj(r) and ad{r) (dotted and wide-dashed 
lines, respectively) that were used by Muchnick & Russek (Ref. 24) are also shown, although their 
ad(r) is not strictly comparable, since they used it in a term with rf rather than f^fa)- 



35 




FIG. 5. Scatterplot of the deviations of the adopted "BMP" He-H 2 PES relative to the 20 203 
ab initio energies, as a function of energy E. Note that for E > 0.13 Eh the vertical scale is 
compressed by a factor of 10; for E > 0.5 Eh, the horizontal scale is compressed likewise. (There 
are 13 397, 4592, and 1991 points in these three parts of the plot, respectively; 223 points lie offscale 
to the right.) Denser bands of points (usually nearly vertical) result from the discrete r values of 
the main ab initio grid. Note that most of the extreme outliers for E > 0.2 Eh arise from the 
conical intersection of the ground state with the first excited state. 
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5 10 5 10 



r (a ) for y = 60° r (a ) for y = 90° 

FIG. 6. Contour plots of the fitted "BMP" He-H2 surface of this paper as a function of H — H 
separation r and H2 — He separation R, at four angles 7 (see Fig. 1). Positive contour levels (solid 
lines) lie at energies E = {0, 0.005, 0.01, 0.02, 0.05, 0.1, 0.2, 0.5, 1, 2, 5, 10, 20, 30, 50, 70, 100, 120, 
140, 160, 180, 200, 300, 400, 500, 750, and 1000 mE h ] relative to separated H+H+He (numbers in 
italics correspond to heavy contour lines — most of these are also labelled in the plots). Negative 
contour levels (dotted lines) lie at energies E = {-0.005, -0.01, -0.025, -0.05, -0.075, -0.1, 
-0.125, -0.15, -0.175, -0.2, -0.25, -0.3, -0.4, -0.5, —1, -2, -5, -10, -20, -40, -60, -80, 
-100, -120, -140, and -160 mE h }. (a) 7 = 0°, (b) 7 = 30°, (c) 7 = 60°, (d) 7 = 90°. Note that 
the small {E ~ —0.3 mEh) "basin" near (r = 10.5 a , R = 0) is spurious: see text. 
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FIG. 7. Comparison of the BMP fit (lines) with ab initio energies (symbols) as a function of 
the H2 — He distance R at 7 = 90° (T-shaped orientation), for five H — H distances r = 1.4, 2.05, 
2.55, 3.2, and 4.0 a . 
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1 2 3 4 5 

Distance R (a ) 



FIG. 8. Potential for the He atom "between" widely-separated H atoms (r = 10.6 a ), for 
displacements R of the He atom at two angles 7 (from the center point). For 7 = 90° (T-shaped: 
heavy lines), as R increases, the He atom is moving out from between the H atoms. For 7 = 0° 
(linear: light lines), as R increases, the He atom approaches one of the H atoms and encounters 
the base of the repulsive "wall." Note that the "modMR" surface is identical to the Muchnick & 
Russek (Ref. 24) surface, except for the long-range (dispersion) terms. 



39 



1.5 

o 

3, 1-4 

| 1-3 

E 1.2 
w 

* 1.1 

cr 

n 
< 

DC 

0.8 



g 0.7 

.52 0.6 
Q 

0.5 



1 

0.9 



.•<•>©• 
.o 



.©©".. 

a ■ 



o . • . . 

• ® ■ ■ .<?'••'•• 

' . ft °. o 
•. » © ■ ■ o© 



O o 1 



o. 



(a) 



dE 1g < 20 mE h 
dE 1g <10mE h 
and < dE 21 / 4 



1.5 2 2.5 3 3.5 4 4.5 5 5.5 6 

Distance r (a ) 



180 
170 
_ 160 



£ 150 

D) 

2- 140 

CD 130 
120 
< 110 

100 



(b) 



dE 1g <40mE h 
dE 1g <20mE h 
dE 1g <10mE h 
and < dE,, / 4 



>•£..•• 



i * * : ■ *■ s » 



1.5 2 2.5 3 3.5 4 4.5 5 5.5 6 

Distance r (a ) 



2.5 - 



(C) 



dE 1g <40mE h 

dE lg <20mE h o 

dE 1g <10mE h □ 

and < dE 21 / 4 ■ 




2.5 



(e) 



dE 1g < 40 mE h 
dE lg < 20 mE, 
dE 1g <10mE h 
and <dE 21 1 A 



1.5 




c 

LJJ 



0.5 1 



1.5 2 2.5 3 3.5 4 4.5 5 5.5 6 6.5 

Distance r (a ) 



1.5 2 2.5 3 3.5 4 4.5 5 5.5 6 6.5 

Distance r (a ) 




1.5 2 2.5 3 3.5 4 4.5 5 5.5 6 

Distance r (a ) 



FIG. 9. The approximate position of the conical intersection between the ground state and 
the first excited state, as a function of the H — H separation r. Plotted points lie close to the 
conical intersection, as they have small differences dE\ g between ground state and first excited 
state energies (20 — 40 mE^ for dots, 10 — 20 mEh for open circles, — 10 mEh for open and filled 
squares — filled squares also have dE\ g < 0.25 where dE<i\ is the energy difference between 

the first and second excited states), (a) The smaller He — H separation Ra, (b) the angle 9 of the 
He atom relative to the closer H atom, (c) the distance R to the center of the H2 molecule, (d) the 
angle 7 relative to the H2 molecule, (e) the energy E relative to separated H + H + He, and (f) the 
error in the fitted BMP surface near the conical intersection. 
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FIG. 10. Four different cuts through (or very near) conical intersection points of the He-H2 
surface, showing ab initio ground state energies (squares) and first excited state energies (dia- 
monds). Smaller symbols indicate points near but not on the cut; their energies have been shifted 
for plotting purposes by the difference in the BMP surface energy corresponding to moving these 
points onto the cut (as expected, this works well for the ground state, but less well for the first 
excited state). Also shown are the fitted BMP surface of this work (solid line) and the previous 
surfaces of Muchnick & Russek (Ref. 24) (dashed line), Wilson, Kapral, & Burns (Ref. 17) (dotted 
line), and Dove & Raynor (Ref. 20) (dot-dashed line). The Schaefer & Kohler surface (Ref. 22) is 
not plotted, as this would require extrapolating it far beyond its range of validity. 
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FIG. 11. Comparison of the fitted BMP surface (solid lines) in the van der Waals well at 
three different orientations with the accurate ab initio energies of Tao (Ref. 31) (circles) and 
the accurate "gen-Tao vdW" energies generated from them (squares); some less-accurate earlier 
ab initio energies (triangles) are also shown, as well as some earlier analytic surfaces (dotted/dashed 
lines: as in Fig. 10). Note the "shifted-logarithmic" energy scale, used to emphasize effects at the 
bottom of the van der Waals well, (a) 7 = 0°, (b) 7 = 45°, (c) 7 = 90°. 
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FIG. 12. Effect of H2-molecule size r on the van der Waals well, at three different orientations 
(note that E R2 = {-83571.4, -171346.4, -174100.0, -167643.3, and -138203.3 /jE h } for r = {0.9, 
1.28, 1.449, 1.618, and 2.0 a }, respectively); the energy scale is the same as in Fig. 11. The fitted 
BMP surface (lines) is compared to the near-H2-equilibrium ab initio energies of Tao (Ref. 31) 
(circles and triangles) and, at more extreme r values, the much-less-accurate generated energies 
described in § II A 4 ("gen-SK vdW": larger crosses and plusses, "gen-MR vdW": smaller ones). 
For the sake of clarity, the "gen- Tao vdW" points are omitted, (a) 7 = 0°, (b) 7 = 45°, (c) 7 = 90°. 
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FIG. 13. Comparison of the fitted BMP surface of this work (and also previous analytic sur- 
faces) to the ab initio energies, for two H2-molecule sizes r. Only for the lowest-energy points does 
the (largely systematic) uncertainty of ~ 0.6 mEh in the ab initio energies approach or exceed 
the size of the symbols, (a) r = 1.4 a , 7 = 0°: near equilibrium H2; the "hole" is visible in the 
Wilson, Kapral, & Burns and Dove & Raynor surfaces (Refs. 17, 20), but the other surfaces are 
quite accurate, (b) r = 4.0 a , 7 = 0°: important for dissociation and high-excitation H2; the 
"hole" in the Wilson, Kapral, & Burns and Dove & Raynor surfaces is at even lower energy, and 
even the Muchnick & Russek surface (Ref. 24) has relatively large fractional errors. The Schaefer 

6 Kohler surface (Ref. 22) is not plotted, as this is too far to extrapolate beyond its range of 
validity, (c) r = 4.0 a , 7 = 45°: similar to (b), but note expanded vertical scale. Smaller symbols 
correspond to points near but not quite on the 7 = 45° cut, shifted as in Fig. 10. (d) r = 4.0 a , 

7 = 90°: similar to (c). 
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